function execute_AVAinv_fixed_dens(ProjectPath, InputPath, OutputPath, XX,YY,ZZ,...
    MinLayerSize,MaxLayerSize, WeightDen, WeightVp, WeightVs, ...
    NbrIt, NbrSim, DenHD, VpHD, VsHD, JointDVp,...
    JointVpVs, KrigAngDen, KrigRangeDen, KrigAngVp, KrigRangeVp, KrigAngVs, KrigRangeVs,...
    VarAngDen, VarRangeDen, VarAngVp, VarRangeVp, VarAngVs, VarRangeVs,...
    ZoneFileName, NbrZones, VariogramType, VarNugget, ...
    SimType, NbrClassesBihist, CovTab, SeisReal, Angles, rangeOfAngles, GatherWavelet, WeightsSelection, ...
    NullDataValue, useLFM, updateLFM, WeightSeismic, ...
    MeanModelFilename)
%% PRE-STACK GLOBAL STOCHASTIC INVERSION WITH MEMORY PERFORMANCE UPDATE. 
% EACH SIM IS SAVED AND NOT STORED IN MEMORY. SAME AS AVO_inv3D.m 
% 
% [FUNCTIONS]
%   runsim_AVO.m        - RUN SIMULATIONS EXTERNALY
%   gslib2grid.m        - CONVERT GSLIB ARRAY INTO GRID
%   fmodel_AVO.m        - FORWARD MODEL FOR AVO INVERSION
%   correlation_AVO.m   - COMPUTE CORRELATION CUBE BTW. ORIGINAL AND SY SEISMIC
%   global_corr_AVO.m   - COMPUTE GLOBAL CORRELATION
%   select_best_AVO.m   - SELECT BEST
%   local_corr_AVO.m    - COMPUTE CORRELATION CUBES BASED ON WEIGHTS
%   grid2gslib.m        - CONVERT FROM GRID INTO GSLIB FORMAT
%
% [MODIFICATIONS]:
% JUL 2019      - Inclusion of fast selection by regression. Cap of CC is
%                 done outside. Handles null data from seismic
% NOV 2018      - Major code update. execution is now done with
%                   execute_AVAinv.
% 21 APR 2014 - MAJOR PERFORMANCE UPDATE
% 03 APR 2014 - PERFORMANCE UPDATE
% 16 DEZ 2013 - LAYERS MAY BE FULLY RANDOMIZED OR NOT - use_layers
% 28 JUN 2013 - PERFORMANCE UPDATE. BEST CUBE FROM PREVIOUS ITERATION IS
%               ONLY KEPT IN MEMORY FOR FIRST ITERATION
% 17 JUN 2013 - SIMULATION ALGORITHM CHANGED AND USER-DEFINED PARAMETERS.
% 12 MAR 2013 - R(0) AND G ARE OPTIONAL CALCULATIONS; AUTOMATIC PATHS
% 17 JAN 2013 - OUTPUT SELECTION
% 23 NOV 2012 - R(0) AND G ARE SAVED WITH THE SYNTHETIC SEISMIC DATA
% 30 OCT 2012 - PROPERTY BOUNDS ARE USER-DEFINED IN THIS SCRIPT
% 03 OCT 2012 - LAYERS ARE COMPUTED AND SAVED PER ITERATION IN MAIN[]
% 28 SEP 2012 - WEIGHTS FOR LOCAL CORRELATION ARE NOW IN MAIN[]
% 11 SEP 2012 - GLOBAL CORRELATION FUNCTION ADDED
% 10 SEP 2012 - LOCAL CORRELATION FUNCTION ADDED. BEST CORRELATION CUBES
%               ARE NOW COMPUTE IN A WEIGHTED VERSION (Vp - NEAR ANGLE,
%               Vs - FAR ANGLE)
% 23 AUG 2012 - layers_sim SAVED
% 21 AUG 2012 - NEW FUNCTIONS runsim.m AND par_file.m
% 20 AUG 2012 - SEED IS SET AS DEFAULT AT THE BEGINNING OF THE INVERSION
%               BEST R(0) AND G ARE WRITTEN
% 07 AUG 2012 - corr_cell IS NO LONGER KEPT IN MEMORY IS SAVED AT EACH
%               ITERATION
% 06 AUG 2012 - SIMULATION ENGINE IS NOW DSS bi-hist WITH SECONDARY IMAGE
% 03 AUG 2012 - SIMULATION GRID; W; y_int; p; LS; interp; - ONLY KEEP 
%               SIMULATIONS FROM CURRENT ITERATION; OUTPUT UPDATED
%               sy_cell SAVED AND CLEARED
% 02 AUG 2012 - NAME OF SIMULATION EMSEMBLE AUTOMATIC
% 01 AUG 2012 - Rzero AND G ARE NO LONGER OUTPUT; LOG FILE UPDATED;
%               SIMULATION ARRAY IS CLEAR AND SAVED
% 31 JUL 2012 - LOG FILE IMPROVED, NOT FINAL YET.
% 18 JUL 2012 - GLOBAL CORRELATION IS NOW USER-DEFINED IN A VARIABLE
% 17 JUL 2012 - MODIFICATIONS FOR PETROBRAS SYNTHETIC
% 12 JUL 2012 - BEST FROM LAST ITERATION IS NOW PART FROM ENSEMBLE SET
% 11 JUL 2012 - global_coor3D.m IS NOW PART OF THE WORKFLOW.
%               QUASI-CORRELATION ADDED FOR FIRST RUNS.
%
% LA - 7 SEP 2012
%
%

% CREATE STRUCTURE FOR SGEMS FILE

if (strcmp(SimType ,'SGEMS') == 1)
    SgemsStrut.z0           = 1;
    SgemsStrut.y0           = 1;
    SgemsStrut.x0           = 1;
    SgemsStrut.xsize        = 1;
    SgemsStrut.ysize        = 1;
    SgemsStrut.zsize        = 1;
    SgemsStrut.nx           = XX(1,1);
    SgemsStrut.ny           = YY(1,1);
    SgemsStrut.nz           = ZZ(1,1);
    SgemsStrut.grid_name    = 'seismic';
    SgemsStrut.property{1}  = 'property';
    SgemsStrut.type_def     = 'Cgrid';
    SgemsStrut.magic_number = 1.561792946e+9;
end

%% INITIALIZE VARIABLES AND ALLOCATE MEMORY
GlobalCorr      = zeros(NbrIt*NbrSim, 1, 'single');
DGridTot        = cell(2,1);
VpGridTot       = cell(2,1);
VsGridTot       = cell(2,1);

%% INITIALIZE RANDOM SEED
rng default

%% Handle hard-data
for i = 1:NbrZones
    
    [~,den_hd_mat]  = import_gslib(DenHD(i,:));
    [~,vp_hd_mat]   = import_gslib(VpHD(i,:));
    [~,vs_hd_mat]   = import_gslib(VsHD(i,:));
  
    % PROPERTY BOUNDS
    BoundsDen(i,:)  = [min(den_hd_mat(:,4)) max(den_hd_mat(:,4))];
    BoundsVp(i,:)   = [min(vp_hd_mat(:,4)) max(vp_hd_mat(:,4))];
    BoundsVs(i,:)   = [min(vs_hd_mat(:,4)) max(vs_hd_mat(:,4))];
    
end
% load mean model if they exist
if(useLFM && ~updateLFM)
    
    [~, DMean]   = import_gslib([ProjectPath,InputPath, MeanModelFilename(1,:)]);
    [~, VpMean]  = import_gslib([ProjectPath,InputPath, MeanModelFilename(2,:)]);
    [~, VsMean]  = import_gslib([ProjectPath,InputPath, MeanModelFilename(3,:)]);
    
    DMean   = gslib2grid(DMean, XX(1,1), YY(1,1), ZZ(1,1));
    VpMean  = gslib2grid(VpMean, XX(1,1), YY(1,1), ZZ(1,1)); 
    VsMean  = gslib2grid(VsMean, XX(1,1), YY(1,1), ZZ(1,1));

end

%%
fileid  = fopen([ProjectPath, OutputPath, 'log.out'], 'w');% INITIALIZE FILE
fprintf(fileid, '%s \n', 'Pre-Stack Global AVA Inversion'); 
fprintf(fileid, '%s \n', '**********************************');
fprintf(fileid, '%s %d \n', 'Number of simulations per iteration', NbrSim); 
fprintf(fileid, '%s %d \n', 'Number of iterations', NbrIt); 
fprintf(fileid, '%s \n', '**********************************');
fprintf(fileid, '%s \n', 'Weights for Density', num2str(WeightDen)); 
fprintf(fileid, '%s \n', 'Weights for Vp', num2str(WeightVp)); 
fprintf(fileid, '%s \n', 'Weights for Vs', num2str(WeightVs)); 
fprintf(fileid, '%s \n', '**********************************');

%% MAIN BODY
for j = 1:NbrIt % CYCLE FOR NO. OF ITERATIONS

    if (useLFM && updateLFM == 1 && j ~= 1)
        DMean   = Dtemp  ./ NbrSim;  
        VpMean  = Vptemp ./ NbrSim;
        VsMean  = Vstemp ./ NbrSim;
        
        grid2gslib(DMean,[ProjectPath, OutputPath,'Dmean_',num2str(j),'.out']);
        grid2gslib(VpMean,[ProjectPath, OutputPath,'VpMean_',num2str(j),'.out']);
        grid2gslib(VsMean,[ProjectPath, OutputPath,'VsMean_',num2str(j),'.out']);
    end
    
    if (useLFM && updateLFM == 1)
        Dtemp   = zeros(XX(1,1), YY(1,1), ZZ(1,1));
        Vptemp  = zeros(XX(1,1), YY(1,1), ZZ(1,1)); 
        Vstemp  = zeros(XX(1,1), YY(1,1), ZZ(1,1));
    end
    
    fprintf(fileid, '%s \n', '**********************************');
    fprintf(fileid, '%s %d \n', 'Iteration number', j); % WRITE FILE WITH NUMBER OF ITERATION
    fprintf(fileid, '%s \n', '**********************************');
    %% CREATE RANDOMIZED LAYERS
    
     layers = make_layers(0,j,OutputPath,MinLayerSize,MaxLayerSize, ZZ(1,1));

    for s = 1:NbrSim % CYCLE FOR NO. OF ITERATIONS       

        %% RUN SIMULATIONS
        [DGridTot{1,1},VpGridTot{1,1},VsGridTot{1,1}] = runsim_AVO_fixed_dens(SimType, ProjectPath, InputPath,...
            OutputPath, j, s, VpHD, VsHD, fileid, ... % eliminated DenHD FT
            BoundsDen, BoundsVp, BoundsVs, XX, YY, ZZ, JointDVp, JointVpVs, KrigAngDen, KrigRangeDen,...
            VarAngDen, VarRangeDen, KrigAngVp, KrigRangeVp, VarAngVp, VarRangeVp,...
            KrigAngVs, KrigRangeVs, VarAngVs, VarRangeVs, VariogramType, VarNugget, ZoneFileName,...
            NbrZones, NbrClassesBihist, CovTab);
        
        if (useLFM && updateLFM == 1 && s == 1)
            
            Dtemp   = DGridTot{1,1};
            Vptemp  = VpGridTot{1,1}; 
            Vstemp  = VsGridTot{1,1};
            
        elseif(useLFM && updateLFM == 1 && s ~= 1)
            
            Dtemp   = Dtemp  + DGridTot{1,1};
            Vptemp  = Vptemp + VpGridTot{1,1}; 
            Vstemp  = Vstemp + VsGridTot{1,1};
            
        end

        tic
        
        syntheticSeis = fmodelFastAVA(2, DGridTot{1,1}, VpGridTot{1,1}, VsGridTot{1,1}, Angles, GatherWavelet);

        twtsim  = toc;
        
        fprintf(fileid, '%s %d \n', 'Forward model - elapsed time (sec)', twtsim); % WRITE NUMBER OF SIM
               
        %% COMPUTE CORRELATIONS
        tic
        if(useLFM && updateLFM == 1 && j ~= 1)
            corr_cell_tot{1,1}      = correlation_AVAmean(1, SeisReal, ...
                syntheticSeis, layers, DMean, ...
                VpMean, VsMean, DGridTot{1,1}, ...
                VpGridTot{1,1}, VsGridTot{1,1}, WeightSeismic);
        elseif(useLFM && ~updateLFM)
            corr_cell_tot{1,1}      = correlation_AVAmean(1, SeisReal, ...
                syntheticSeis, layers, DMean, ...
                VpMean, VsMean, DGridTot{1,1}, ...
                VpGridTot{1,1}, VsGridTot{1,1}, WeightSeismic);
        else
            corr_cell_tot{1,1}      = correlation_AVA(0, SeisReal, syntheticSeis, layers);
        end
        twtsim = toc;
        fprintf(fileid, '%s %d \n', 'Correlations - elapsed time (sec)', twtsim); % WRITE NUMBER OF SIM
        
        %% COMPUTE GLOBAL CORRELATION
        tic
        
% FT original        GlobalCorr = global_corr_AVA(SeisReal, syntheticSeis, GlobalCorr, j, s, fileid, NbrSim, NullDataValue);
            GlobalCorr = global_corr_AVA(SeisReal, syntheticSeis, GlobalCorr, j, s, fileid, NbrSim, NullDataValue, rangeOfAngles);   
            
        twtsim = toc;
        fprintf(fileid, '%s %d \n', 'Global Correlation - elapsed time (sec)', twtsim); % WRITE NUMBER OF SIM
        
        clear sy_seis
                
        % FIRST SIMULATION BOTH ARRAYS (SIM {1,1} AND BEST {2,1} ARE THE SAME
        if s==1 
            tic
            
            corr_cell_tot{2,1}  = corr_cell_tot{1,1}; 
            DGridTot{2,1}       = DGridTot{1,1}; 
            VpGridTot{2,1}      = VpGridTot{1,1}; 
            VsGridTot{2,1}      = VsGridTot{1,1};
            
            twtsim=toc;
            fprintf(fileid, '%s %d \n', 'Allocating arrays - elapsed time (sec)', twtsim); % WRITE NUMBER OF SIM
        else
            %% SELECT BEST PER ITERATION - FROM THE REALIZATION SET
            % ITERATION 1++ THERE IS BEST IN {3,1}

            fprintf(fileid, '%s %d \n', 'Select Best - Simulation',s);
            tic
            
%             [DGridTot,VpGridTot,VsGridTot,corr_cell_tot]=...
%                 selectBestAVAangles(corr_cell_tot,DGridTot,VpGridTot,VsGridTot);
%             [DGridTot,VpGridTot,VsGridTot,corr_cell_tot]=...
%                 select_best_AVA_geom(corr_cell_tot, Angles, layers, DGridTot, VpGridTot, VsGridTot);  
            [DGridTot,VpGridTot,VsGridTot,corr_cell_tot] =...
                select_best_AVA_regression(Angles, rangeOfAngles, WeightsSelection,corr_cell_tot, layers, DGridTot, VpGridTot, VsGridTot);
            twtsim=toc;
            fprintf(fileid, '%s %d \n', 'Best fit selection - elapsed time (sec)', twtsim); % WRITE NUMBER OF SIM
        end
    end
    tic
    %% COMPUTE LOCAL CORRELATIONS
    tic
    corr_cell_tot{2,1} = capCC_AVA(j, corr_cell_tot{2,1});
    
    %     localCorrAVAangles(SimType, ProjectPath, OutputPath, j, corr_cell_tot, SgemsStrut)    
    local_corr_AVA(SimType, OutputPath, j, layers, corr_cell_tot, XX(1,1), YY(1,1), ZZ(1,1), WeightDen, WeightVp, WeightVs, SgemsStrut, NullDataValue)    
    
    twtsim=toc;
    fprintf(fileid, '%s %d \n', 'Computing and Writing BCDen, BCVp, BCVs- elapsed time (sec)', twtsim); % WRITE NUMBER OF SIM
    
    %% CREATE AND WRITES GLOBAL BEST = BEST IT 1
    
    % assign nulldatavalues to NaN
    DGridTot{2,1}(find(isnan(DGridTot{2,1}) == 1))      = NullDataValue;
    VpGridTot{2,1}(find(isnan(VpGridTot{2,1}) == 1))    = NullDataValue;
    VsGridTot{2,1}(find(isnan(VsGridTot{2,1}) == 1))    = NullDataValue;
    
    if (strcmp(SimType ,'SGEMS') == 1)
        SgemsStrut.grid_name        = 'Density';
        SgemsStrut.data             = DGridTot{2,1}(:);
        SgemsStrut.property{1,1}    = 'Density';
        sgems_write([OutputPath,'best_D_',num2str(j),'.sgems'], SgemsStrut);
        
        SgemsStrut.grid_name        = 'Vp';
        SgemsStrut.data             = VpGridTot{2,1}(:);
        SgemsStrut.property{1,1}    = 'Vp';
        sgems_write([OutputPath,'best_vp_',num2str(j),'.sgems'], SgemsStrut);
        
        SgemsStrut.grid_name        = 'Vs';
        SgemsStrut.data             = VsGridTot{2,1}(:);
        SgemsStrut.property{1,1}    = 'Vs';
        sgems_write([OutputPath,'best_vs_',num2str(j),'.sgems'], SgemsStrut);

    else
        % SAVE BEST CUBES FROM CURRENT ITERARION
        grid2gslib(DGridTot{2,1}, [OutputPath,'best_D_',num2str(j),'.out']);
        grid2gslib(VpGridTot{2,1},[OutputPath,'best_vp_',num2str(j),'.out']);
        grid2gslib(VsGridTot{2,1},[OutputPath,'best_vs_',num2str(j),'.out']);
    end
    
    % ASSIGN FINAL BEST PER ITERATION TO {3,1}
    twtsim=toc;
    fprintf(fileid, '%s %d \n', 'Writing best models - elapsed time (sec)', twtsim); % WRITE NUMBER OF SIM
end

%% Save global CC evolution
save([ProjectPath, OutputPath,  'globalCorr.mat'], 'GlobalCorr');

fclose(fileid); % CLOSE LOG FILE
end